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Phototactic Clustering of Swimming Micro-organisms in a Turbulent Velocity Field 
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We study the distribution of swimming micro-organisms advected by a model turbulent flow and 
attracted towards a localised light source through phototaxis. It is shown that particles aggregate 
along a dynamical attractor with fractal measure whose dimension depends on the strength of the 
phototaxis. Using an effective diffusion approximation for the flow we derive an analytic expression 
for the phototactic gain (increase in light exposure over the aggregate) and by extension an accurate 
prediction for the fractal dimension based on the properties of the advection dynamics and the 
statistics of the attracting field. This shows that the fractal characteristics of the aggregate are 
determined by the non-dimensional ratio of the kinetic energy of swimming to that of the turbulent 
flow. 

PACS numbers: Valid PACS appear here 



The inhomogeneous distribution of chaotically ad- 
vected particles has received much recent attention in 
a variety of contexts. It has been shown in the case of 
passive particles in compressible surface flows lEIS Hi 
and inertial particles in incompressible flow 0, HQ that 
clusters form along a dynamical attractor of fractal mea- 
sure. This intermittent distribution has important con- 
sequences for physical process such as rain formation Q 
and chemical or biological reaction processes [1,01 . In the 
context of biological processes microscale patchiness in 
plankton populations has been observed in several stud- 
ies [HI HI 13, HI and is seen to play an important role 
in plankton population dynamics and therefore oceanic 
ecology. 

In this letter we consider swimming particles advected 
by an incompressible flow with an effective compressibil- 
ity induced by particle motility in the presence of an at- 
tracting light source. Particles are able to detect the 
gradient of the illumination field and move in the direc- 
tion of increasing light intensity with a swimming veloc- 
ity assumed to be proportional to the detected gradient. 
Motion is therefore governed by the following equation 



r = v/(r,t) + xV$(r) 



(1) 



where the parameter x is the phototactic coefficient 
which defines the strength of the particles reaction to 
the gradient of the illumination field ^(r). For our simu- 
lations we assume a Gaussian light distribution although 
results are not dependent on the exact functional form of 
$(r). 

As a carrier flow we use a two-dimensional synthetic 
turbulence model which approximates a turbulent flow 
lj,ll5|. It is Gaussian, isotropic and statistically homo- 
geneous with a prescribed energy spectrum of the form 



originally proposed by Kraichnan 16]. The flow is gener- 
ated by using an Ornstein-Uhlenbeck process defined as 
the solution to the stochastic partial differential equation 



(3) 
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where ip is the stream function of the velocity field, v is a 
viscosity parameter that sets the time correlation of the 
flow, £ represents the magnitude of the stochastic forcing 
and W is a Wiener process. Following [lH we introduce 
the lengthscale L and timescale T such that T = Lnjv~j\ 
and rescale so that the flow is periodic on a square domain 
with unit length. The lengthscale of the highest energy 
modes is defined by 27r/fco and we set this value to be two 
thirds of the box size. The spectrum is normalised which 
leads to a mean squared velocity of the flow of 0.5. For 
the illumination field we use a Gaussian function with an 
amplitude of unity and a width of a quarter of the box 
length, $(r) = exp [-8 (x 2 + y 2 ) / L 2 ] . 

In an ergodic incompressible flow with no transport 
barriers an ensemble of passively advected particles are 
quickly dispersed producing a spatially uniform density 
distribution. Conversely phototactic particles swimming 
in the direction of a fixed illumination gradient in the ab- 
sence of advection would concentrate in regions of max- 
imum light intensity. When combined these two effects 
result in a statistically stationary steady state particle 
distribution which is highly non-uniform in space. Par- 
ticle distributions after the relaxation of transients are 
shown in this regime in Fig. [1] Particles clearly do not 
converge to a fixed point as would be the case for strong 
phototactic attraction but instead aggregate along a dy- 
namically evolving fractal attractor and continue to visit 
all regions of the physical space. 

Despite this, the average illumination received by the 
ensemble of particles (or equivalently along the trajec- 
tory of a single particle) increases with \ an d is higher 
than the spatial average illumination. In Fig. [2] we show 
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the time-averaged increase in light exposure of an en- 
semble of phototactic particles as compared to passively 
advected non-motile particles as a function of the photo- 
tactic coefficient. 

As x is increased the attractor becomes more singular 
with most particles concentrated in thin filaments sepa- 
rated by empty regions. To quantify the characteristics 
of the attractor we numerically calculate the information 
dimension defined as 



17] 
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(4) 



where the domain is covered by a rectangular grid of box 
size e and pi is the fraction of particles in a given box i. 
Results for various values of \ are shown in Fig. [51 

In order to investigate the relationship between the 
phototactic coefficient x an d both the average illumina- 
tion gain of the aggregate and the fractal dimension of 
the attractor we approximate the turbulent advection as 
a purely diffusive process. 

Although the instantaneous dynamical attractor 
changes in an irregular fashion following the turbulent 
flow, by averaging over time we can obtain a picture of 
the denser regions of the physical space which will then 
allow us to quantify various statistical properties of the 
system. To this end we define p as the time averaged den- 
sity of particles at a given point in space and the notation 
(..)p to mean an average quantity defined over p. With 
these definitions and an effective diffusion approximation 
for the flow we model the system by the equation 



dp 

dt 



D f V 2 p- 



V • (pxV$) 



(5) 



where D j is the effective diffusivity of the flow defined by 
the asymptotic dispersion rate of advected passive parti- 
cles in an unbounded domain as (d 2 ) ~ 4Z?/i. The first 
term on the r.h.s. of Eqn. [5] represents the diffusive ef- 
fect of the turbulent flow which smooths out any inhomo- 
geneities in the distribution and if dominant would lead 
to a uniform particle density. The second term counter- 
acts this effect by collecting particles in high light inten- 
sity regions leading to a non-uniform distribution. Eqn. [5] 
has a steady state solution of the form 



P(r) = i exp 



(6) 



which is analogous to a Boltzmann distribution of non- 
interacting particles in an external field, where Z is a 
normalisation factor akin to a partition function which 
allows us to use p as a probability density 



exp 



Di 



$(r) 




dV. 



(7) 



FIG. 1: Distribution of 500,000 phototactic particles after 
relaxation of transients. Above: \ = 0.15. Below: \ = 0.05. 
Color version: Light distribution is shown in yellow 



To compare this result with the numerical experiments 
we first note that as a result of the radial symmetry of 
the illumination field the particle distribution is also ra- 
dially symmetric. This allows us to fully describe the 
two dimensional distribution by a one dimensional radial 
profile. In Fig. [3] the result predicted from Eqns.[6]and 
[7] is compared to a time averaged, normalised histogram 
plot of radial location of particles for x = 0.05. We also 
include in the inset of Fig. [3] a scatter plot of numerical 
results against expected values for a range of x which 
illustrates the excellent agreement of numerics with the 
effective diffusion approximation. 

Using the distribution given above we can evaluate the 
average light intensity experienced by particles with a 
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FIG. 2: Above: Gain in light exposure achieved by photo- 
taxis. Squares - numerical results, solid line - theoretical 
approximation, dashed line - first order theoretical approx- 
imation. Below: Information dimension. Squares - values 
calculated from distribution by Eqn. 2] dashed line - results 
from numerically calculated Lyapunov exponents, solid line - 
theoretical model. 
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FIG. 3: Radial plot of probability density function (solid line) 
and numerical results (points) for average particle distribu- 
tion, x — 0.05. Inset: Scatter plot of predicted density (p-r) 
to numerically observed density (pjv) for range of x from 0.025 
to 0.15. Points lying on the diagonal represent exact match 
between effective diffusion model and numerics. 



given phototactic response parameter \, as 



\ / P 7 



exp 



which can be rewritten as 



7T* (r) 
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dx 



(8) 



(9) 



Then as Z is an expectation value over the unit area do- 
main the logarithmic term can be treated as a cumulant 
generating function of the light distribution <£>. This leads 
to 



71=0 
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(10) 



where the cumulants of the illumination field. For 

weakly phototactic particles, i.e. x/Df <C 1, the higher 
order terms can be neglected leaving 



(ii) 



This explains the linear relationship seen in Fig. [2] for 
small values of x an d shows that the gain in illumina- 
tion exposure due to phototaxis is proportional to the 
variance of the illumination field. The truncated series is 
plotted in Fig. [2] alongside results containing higher order 
terms required for convergence. 

Interestingly, the time-averaged distribution defined by 
Eqn. [6] can also be used to calculate the dimension of the 
dynamical attractor. The Kaplan- Yorke conjecture [l8| 
relates the information dimension of an attractor to the 
average stretching and contraction rates in the phase 
space by the formula (for two dimensional phase space) 



Di = 1 



Ao_ 

|A 2 | 



(12) 



where Ai, A2 are the Lyapunov exponents of the system. 
The sum of the Lyapunov exponents defines a dissipa- 
tion rate a, i.e. the average rate at which volumes in the 
phase space contract. In a conservative system as is the 
case of particles passively advected by an incompressible 
flow a is zero and therefore in the two dimensional model 
Ai = — A 2 and the information dimension is two. How- 
ever if phototaxis is present particles no longer follow the 
streamlines of the flow and there is a net contraction in 
the system which reduces the dimension of the attractor. 

The Lyapunov exponents can be determined numeri- 
cally by a standard method proposed by Benettin [T3| 
and the values found in this way show the dependence of 
Ai, A2 and a on the phototactic coefficient (Fig.|4|). The 
information dimension based on these quantities can now 
be found by using the Kaplan- Yorke conjecture and this 
agrees with the values calculated directly from the par- 
ticle distributions as shown in Fig. [2] 

Returning to the theoretical model the dissipation rate 
a can be obtained from the divergence of the velocity field 
of the phototactic particles. Due to the incompressibility 
of the carrier flow this gives 



( X V 2 $)_= / p X V 2 <P(r)dV 



(13) 



where by ergodicity averaging over a single trajectory 
over long times is taken to be equivalent to the weighted 
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FIG. 4: Lyapunov exponents and dissipation rate. Symbols 
represent numerical values found by the Benettin method, 
lines values calculated from the theoretical approximation de- 
fined by Egns. [131 and [T4l : Squares/solid line - An,, crosses/ 
heavy dashed - a, asterisks/fine dashed - A2. 



average over the distribution p previously defined in our 
continuum description. 

From the isotropic properties of the flow and zero cor- 
relation between fluid dynamics and illumination field 
we expect the contraction of phase space elements due 
to phototaxis to have an equal effect on both Lyapunov 
exponents. This is consistent with numerical results as 
shown in Fig. 0J Defining Ao as the positive exponent 
of the carrier flow (i.e. in the case of passive particles 
Ai = Ao, A2 = — Ao) the exponents for the phototactic 
particles are then given by 

Ai = A + i ( X V 2 $) _ A 2 = -Ao + \ (xV 2 $)_ (14) 

For the Gaussian light distribution the integral in Eqn. 1131 
can be evaluated numerically and values for Ai, A2 and 
a are plotted in Fig. [4] alongside those obtained from the 
Benettin method. Combining Eqns. [T^l [32] and Q3] gives 
an expression for the fractal dimension of the distribution 
in terms of <£>, \ and the properties of the carrier flow field 
Aq and Df. 



1 



1 



2 4A Z 



exp 



X V 2 $dV 



(15) 



To illustrate their agreement with the numerical results 
solutions to this equation are included in Fig. O 

An alternative formulation can be obtained from 
Eqn. [13] by the use of integration by parts 



-g- (xV$) 2 dV (16) 



where the boundary term vanishes for most boundary 
conditions (i.e. smooth periodic or no-flux illumination 



field) [20|]. The second integral demonstrates that a is 
always negative and is proportional to the mean square 
swimming velocity of the phototactic particles. Thus the 
information dimension can be written as 



£>i = 2 



(VI) 
2X D f 



(17) 



In a fully chaotic flow the effective diffusivity can be es- 
timated as Df ~ UL and the Lyapunov exponent is an 
inverse characteristic timescale of the flow Ao ~ U/L. 
Using these relations Eqn. [T71 shows that the information 
dimension is a function of the non-dimensional ratio of 
the average kinetic energy of swimming measured relative 
to that of the carrier flow. 

We have therefore seen how the effects of phototaxis 
on the spatial distribution and increase in light exposure 
of an ensemble of micro-organisms can be explained by 
an analysis of the attracting field and the characteristics 
of the carrier flow. Specifically we have related the di- 
mensionality of the distribution to the measurable kinetic 
energy of the swimming organisms. 

The combined effect of phototaxis and turbulent 
flow leads to highly non-uniform distribution of micro- 
organisms and this provides a plausible mechanism for 
the generation of microscale patchiness observed in ex- 
periments. Our analytical results for the key character- 
istics of systems of phototactic particles may be used to 
make testable quantitative predictions or to extract infor- 
mation from measurements (e.g. based on the relation- 
ship between swimming kinetic energy and dimension). 
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